function [elements,processInfo] = CalculateMaterialResponse...
    (elements,model,processInfo)
ne = model.countElements();
dim = 2;
for i=1:ne
    e = model.getElement(i-1);
    f = e.getFace();
    ue = f.getUHatLocal();
    nn = f.countModes();
    ng = e.getNumGaussPoint();
    C = elements(i).constitutiveMatrix;
    alpha = elements(i).alpha;
    
    for j=1:ng
        xi = e.getGaussPoint(j-1);
        w = e.getGaussWeight(j-1);
        jac = f.jacXAt(xi);
        detJac = det(jac);
        %compute strain
        B = zeros(3,dim*nn);
        B = e.computeB(f,B,xi);
        G = computeEnhancedStrainOperator(e,xi);
        epsilon = B*ue + G*alpha;
        %compute stress
        sigma = C*epsilon;
        %compute h
        if j==1
            h = w*G'*sigma*abs(detJac);
        else
            h = h + w*G'*sigma*abs(detJac);
        end
%         elements(i).h = h;
        elements(i).stress(:,j) = sigma;
    end
    %update h
    elements(i).h = h;
    %update alpha
    H = elements(i).H;
    Gamma = elements(i).Gamma;
    deltau = ue - elements(i).ue_old;
    alpha = alpha - H\(Gamma*deltau - elements(i).h);
    elements(i).alpha = alpha;
    %update ue_old
    elements(i).ue_old = ue;
end
